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We present a new algorithm for evolving orbiting black-hole binaries that does not require excision 
or a corotating shift. Our algorithm is based on a novel technique to handle the singular puncture 
conformal factor. This system, based on the BSSN formulation of Einstein's equations, when used 
with a 'pre-coUapsed' initial lapse, is non-singular at the start of the evolution, and remains non- 
singular and stable provided that a good choice is made for the gauge. As a test case, we use 
this technique to fully evolve orbiting black-hole binaries from near the Innermost Stable Circular 
Orbit (ISCO) regime. We show fourth order convergence of waveforms and compute the radiated 
gravitational energy and angular momentum from the plunge. These results are in good agreement 
with those predicted by the Lazarus approach. 
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One of the most significant goals of numerical relativity 
is to compute accurate gravitational waveforms from as- 
trophysically realistic simulations of merging black-hole 
binaries. The expectation of very strong gravitational 
wave emission from the merger of two black holes, and 
some of the newest astrophysical observations jfrom su- 
permassive galactic nuclei just about to merge 111 to stel- 
lar size black-hole binaries, make these systems one of the 
most extraordinary astrophysical objects under study to- 
day. Binary black hole mergers are expected not only to 
provide information about the history and formation of 
the binary system but also to provide important precise 
tests of strong-field, highly dynamical relativity. 

Motivated by the forthcoming observations of ground- 
based gravitational wave detectors, such as LIGO 2] , and 
by the next generation of space-based detectors, such as 
LISA [3, the numerical relativity community has dedi- 
cated a great deal of effort to solving the binary-black- 
hole problem over the past few decades. After the 'Bi- 
nary Black Hole Grand Challenge' Q several new ap- 
proaches have been pursued in the attempt to produce 
stable three-dimensional (3D) numerical codes capable of 
evolving the full Einstein field equations in the absence of 
any symmetry. This includes the introduction of new for- 
mulations of these equations and the development of nu- 
merical techniques for accurate evolutions of black-hole 
binaries, such as higher order finite differencing, spectral 
methods, and adaptive mesh refinement (AMR) (see Ref. 
[3 and references therein). 

The calculation of the gravitational radiation emitted 
from plunging black-hole binaries was pioneered through 
the use of the Lazarus approach, which bridges numerical 
relativity and perturbative techniques to extract approx- 
imate gravitational waveforms [a> S ■ More recently 
important progress has been made toward evolving orbit- 
ing binary-black-hole spacetimes with the use of stable 
full 3D numerical relativity codes using corotating gauge 



conditions and singularity excision 




Here we present a novel technique for evolving orbit- 
ing black holes based on puncture data. This technique 
does not require a corotating shift or singularity exci- 
sion. Most importantly, we can produce accurate com- 
plete waveforms from merging black-hole binaries. 

In a previous paper we presented techniques for suc- 
cessfully performing numerical relativity simulations of 
black-hole binaries with fourth-order accuracy Our 
simulations are based on a new coding framework, LazEv, 
which is built on top of the Cactus Computational 
Toolkit 'is*!; this currently supports higher order finite 
differencing for the BSSN formulation of Einstein's equa- 
tions 0, EE 01 1 but is designed to be readily applicable 
to a broad class of formulations. Highly accurate evo- 
lutions can be achieved using our unigrid higher-order 
finite-difference code along with a non-uniform coordi- 
nate system, such as 'Fisheye' 8j, that concentrates grid- 
points in the central region containing the black holes. 

In the puncture ap proa ch |l7l | the metric on the ini- 
tial slice is given by [l^ 7ab = (V'bl + uY^ab, where 
i'BL = 1+X]r=i '^'Z Brill-Lindquist conformal 
factor, mi is the mass parameter of puncture i, is the 
coordinate distance to puncture i, and u is finite on the 
punctures. If the puncture positions are fixed throughout 
the evolution, then the singular behavior in the metric is 
contained 'm ipBL and can thus be treated analytically. 
However, holding the puncture fixed throughout the evo- 
lution leads to significant coordinate distortions that tend 
to kill the run before a common horizon forms. 

We propose a method for evolving puncture type data 
without fixing the puncture positions during the evo- 
lution. Our method is based on the BSSN formula- 
tion. In the BSSN system one evolves a conformal 
metric 7ah = cxp(— 4(/))7ab which has unit determinant, 
K — K"^^, the conformal trace-free extrinsic curvature 
Aab = <2-Kp{—A(l)){Kab — JabK/3), the conformal expo- 
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nent 4>^ a-nd = —dj^^^. In order to regularize the sys- 
tem near the puncture, we replace the BSSN exponent (p 
(which has an O(lnr) singularity at the puncture) with 
a new variable x — exp(— 4(/)) that is on the punc- 
ture. Additionally we modify the standard 1-l-log lapse 
and Gamma-driver shift gauge conditions 0|. Our sys- 
tem is explicitly finite on the initial slice. The evolution 
equations for this system are 



dtx = \x{o^K~dan + P'da, 

doAij = x{-DiD.ja + aRij) 



TF 



doK = ^DW,a + a { A, A' , ^ 



-K- 



(1) 
(2) 

(3) 
(4) 



dtT 



2a ( r',fei^'= + 6i*^5,0 - -f'djK ] , (5) 



where do — dt — Cp, TF indicates that only the trace-free 
part of the tensor is used, Rij — Rij + R^ij is given by 



Ri 



-2D,DjC^ - 2j,jD^Dk(f> + AD,c^D,(f> - 
A^,,b''cpbk<t>, (6) 



(7) 



Di is the covariant derivative with respect to 7ij, and 
Di is the covariant derivative with respect to 7ij. is 
replaced by —dj-y^^ in Eqs. - tZI wherever it is not 
differentiated. Note that Eqs. lj2Jl and (O give the dt 
derivatives of x a-nd F* rather than the do derivatives. 
The Lie derivatives of the non-tensorial quantities (7^ , 
and Aij) are given by 

7,fe5,/3'^- - ^7y■9fe/3^ (8) 
CpA^j = f3''dkA,j+A,^djl3^ + 

Ajkd^p"" - '^A^.dkP". (9) 

Note that dicj) = -l/(4x)9iX and = \{-dzjx/x + 
dixdjx/x^) are singular on the puncture; as a result, we 
identify several potentially singular terms in Eqs. O-©. 
In Eq. © the term D,D ja can be expressed in terms 
of di(f>dja + non-singular terms. However this term is 
multiplied by x and the product is on the puncture. 
Additionally, the didjcf) terms (from Rij) are multiplied 



by ax, and are thus on the puncture. If the lapse 
a ~ on the puncture (as is our choice) then Aab is 
on the puncture (provided that it is on the initial 
slice). However, in Eq. J^J we find the singular term 



aA^^di(j). With our choice of initial data, A^^ 



and 



thus F* is C° on the puncture. We can then choose an 
initial lapse a — V'bl which is O(r^) on the puncture. 
With this choice of lapse, F' evolves to a function that 
is on the puncture. However, F' will only remain 
well behaved if the lapse condition maintains this O(r^) 
behavior near the puncture. We found that the following 
gauge conditions produced smooth waveforms: 

doa = -2aK, (10) 
dtl3'' = B", dtB" = 3/Adtf'' -vB". (11) 

In black-hole binary systems the features that need 
to be resolved range in scale from a fraction of M near 
the horizons to over lOOM in the wave zone when black 
holes are in a slow inspiral motion. This difference in 
scales makes simple unigrid evolutions extremely ineffi- 
cient. We mitigate this problem by introducing a 'multi- 
ple transition' Fisheye transformation. This Fisheye co- 
ordinate is a natural extension of the 'transition Fisheye' 
coordinate 0, 0| (where the effective resolution changes 
from some inner resolution h to an outer resolution a h 
inside a region of a given width, where a is a user specified 
parameter), but with multiple transition regions allowing 
for fine tuning of the resolution in intermediate regions. 
This new Fisheye coordinate simulates fixed-mesh refine- 
ment. 

For puncture data the estimated ISCO ,2^| is charac- 
terized by the parameters 



L/M = 



4.9, P/M ^ 0.335, 
= 0.77, Mn = 0.178, 



Y/M = ±1.1515, 
m = 0.45Af, (12) 



where m is the mass of each single black hole, M is the 
total ADM mass of the binary system, L is the proper 
distance between the apparent horizons, P is the mag- 
nitude of the linear momenta (equal but opposite and 
perpendicular to the line connecting the holes), J is the 
total angular momentum, and (0, Y, 0) is the coordinate 
location of the punctures. We use the Brandt-Briigmann 
approach along with the BAM_Elliptic [isl [l8^ Cactus 
thorn to solve for these initial data. 

We use TT-rotational symmetry about the polar axis 
(Pi-symmetry) and reflection symmetry across the or- 
bital plane to reduce the computational domain to one 
quadrant. We choose Fisheye parameters that produce 
an inner resolution h, an intermediate resolution of 6/1, 
and an outer resolution of 25ft- (to push the boundaries 
out to 60M). 

At the puncture x = and Eq. |(5J| implies that the 
puncture position obeys dtXpunct = -ffixpunct)- We use 
this to track the puncture positions throughout the evo- 
lution. Figure^shows the trajectory of the punctures for 
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FIG. 1: The trajectories of the punctures (the circles cor- 
respond to positions at t=0, 2.5M, 5M, etc). The common 
horizon forms just after the puncture complete a half orbit. 
The punctures continue to orbit throughout the evolution. 
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FIG. 2: QCO waveform. The real and imaginary parts of the 
{£ — 2,m — 2) mode of ^pi calculated at r = lOM. 



an h = M/24 run. A common horizon forms just after 
the punctures complete half of an orbit. The punctures 
continue to orbit throughout the evolution. 

We use the Zorro thorn ^] to calculate ipi and 
decompose it into {i,m) modes. Figure |21 shows the real 
and imaginary parts of the (£ — 2,m — 0) mode at r = 
lOM. We stopped this run at ^ 85M because of limited 
computational resources. Figure shows the real part of 
the {£ — 2,m = 2) mode of -04 at r = 5M (we choose this 
small observer radius to delay outer boundary effects) 
for resolutions of h = Af/16, M/24, M/36, as well as a 
convergence plot of these data. The convergence plot 
shows that the waveforms are fourth-order convergent. 
The radiated energy (as measured from the M/24 run) is 
2.8% ±0.2% in excellent agreement with the final horizon 
mass (see below) , and the radiated angular momentum is 

= -O.llM^iO.OlM^ (computed using Eqs. (22)-(24) 
of Ref. 113). 

We use Jonathan Thornburg's AHFinderDirect 
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FIG. 3: QCO waveforms. The top plot shows the real part of 
the (^ = 2, m = 2) mode at r = 5M for resolutions oi h = 
M/16, M/24, M/36. The bottom plot shows the differences 
between waveforms for h — M/16 and h = M/24: as well 
as the difference between waveforms for h = M/24: and h = 
M/36. The latter difference has been rescaled by (3/2)* to 
demonstrate fourth-order convergence. 



thorn 22] to find apparent horizons. We first detect a 
common apparent horizon at t — 19M. The common 
horizon has an irreducible mass of 0.9055M and the 
ratio of polar to equatorial circumferences asymptotes 
to 0.898 ± 0.002. One can show analytically % that for 
a Kerr black hole, the ratio of the polar and equatorial 
horizon circumferences Cr = Cp/Ce is given by 



i + Vi 



■E 



(13) 



where a — a/Mji and E(x) is the complete elliptic in- 
tegral of the second kind. In the case of a perturbed 
black hole produced by a merger, this ratio shows quasi- 
normal ringing behavior before damping to the expected 
Kerr value. The horizon mass is related to the spin and 

irreducible mass by Mu = (Mii.r/a)^2(l — Vl — a^). 
Hence, the irreducible horizon mass and circumference 
ratio that we measure correspond to a spin of a = 
0.683 ± 0.006 and a horizon mass of 0.973M ± 0.002M. 
The horizon mass reduction is in excellent agreement 
with the calculated radiated energy of 0.028M±0.002A/. 

Table m summarizes the main results of our full numer- 
ical evolution of binary black holes from the ISCO down 
to the final Kerr black hole remnant. Waveforms, as de- 
scribed in terms of the Weyl scalar "04, are dominated 
by the modes £ — 2,m = ±2 and show a strong circular 
polarization as seen along the axis of orbital symmetry. 
These results are in good agreement with the results cal- 
culated by the Lazarus approach 

The Hamiltonian constraint violation along the y axis 
(when the punctures are no longer on the axis) shows cu- 
bic convergence in the former positions of the punctures 
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TABLE I: Results of the evolution. 



Method 


Erad/M 


Jrad/M'^ Merger 


time/M 


a/A'/7i 


This paper 


2.8 ± 0.2 


14 ± 1% TCAH ? 


a 19 


0.683 ± 0.006 


Lazarus" 


2.5 ±0.2 


13 ± 2% Txran ^ 


s 10 


0.70 ± 0.02 



"Errors quoted in the Lazarus runs arc only those from the dif- 
ferences among transition times; hence they only represent a lower 
bound to the total errors. 



and quadratic convergence outside. This quadratic con- 
vergence is most likely due to our second-order-accurate 
initial data. Further study is needed to see if constraint 
violations on the puncture converge to zero. 

Note that in this paper, we used a convenient numeri- 
cal tetrad to calculate the Weyl scalar i/j^. This waveform 
extraction procedure is valid when the far-field spacetime 
approaches that of a perturbed Schwarzschild black hole. 
In future simulations we plan to use a more robust wave 
extraction method recently implemented in Ref. |23l |. 
This method is based on the use of a quasi-Kinnersley 
tetrad to compute -04 and should allow us to extract 
waveforms closer to the hole, where outer boundary ef- 
fects are less important. 

The full nonlinear numerical technique just described 
has shown long-term stability, at least for the merger 
of binary black hole cases we have had the opportu- 
nity to analyze so far. We have also been able to prove 
fourth-order convergence up to relatively high resolutions 
such as h = M/36. This opens up the possibility of 
studying even more interesting astrophysical scenarios 
such as black-hole binaries starting from larger separa- 
tions, which would undergo several orbits before the fi- 
nal plunge. Additionally, we plan evolutions based on 
the more astrophysically relevant Thin-Sandwich and 
Post-Newtonian initial data sets. We will also con- 
sider unequal-mass black-hole binaries and compute their 
gravitational kick in order to evaluate its astrophysical 
consequences [2^ . Finally, we plan to examine merg- 
ers of highly spinning black-hole binaries and study the 
possible 'hang-up' of the binary until the excess angular 
momentum is radiated. 

This work was recently presented at the "Numerical 
Relativity 2005: Compact Binaries" i27|. See [13 for a 
similar approach with comparable results. 
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at Brownsville (NAG5-13396) and the NSF for financial 
support from grants PHY-0140326 and PHY-0354867. 
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